home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / zlahqr.f < prev    next >
Text File  |  1996-07-19  |  12KB  |  381 lines
  1.       SUBROUTINE ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
  2.      $                   IHIZ, Z, LDZ, INFO )
  3. *
  4. *  -- LAPACK auxiliary routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     September 30, 1994
  8. *
  9. *     .. Scalar Arguments ..
  10.       LOGICAL            WANTT, WANTZ
  11.       INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
  12. *     ..
  13. *     .. Array Arguments ..
  14.       COMPLEX*16         H( LDH, * ), W( * ), Z( LDZ, * )
  15. *     ..
  16. *
  17. *  Purpose
  18. *  =======
  19. *
  20. *  ZLAHQR is an auxiliary routine called by ZHSEQR to update the
  21. *  eigenvalues and Schur decomposition already computed by ZHSEQR, by
  22. *  dealing with the Hessenberg submatrix in rows and columns ILO to IHI.
  23. *
  24. *  Arguments
  25. *  =========
  26. *
  27. *  WANTT   (input) LOGICAL
  28. *          = .TRUE. : the full Schur form T is required;
  29. *          = .FALSE.: only eigenvalues are required.
  30. *
  31. *  WANTZ   (input) LOGICAL
  32. *          = .TRUE. : the matrix of Schur vectors Z is required;
  33. *          = .FALSE.: Schur vectors are not required.
  34. *
  35. *  N       (input) INTEGER
  36. *          The order of the matrix H.  N >= 0.
  37. *
  38. *  ILO     (input) INTEGER
  39. *  IHI     (input) INTEGER
  40. *          It is assumed that H is already upper triangular in rows and
  41. *          columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1).
  42. *          ZLAHQR works primarily with the Hessenberg submatrix in rows
  43. *          and columns ILO to IHI, but applies transformations to all of
  44. *          H if WANTT is .TRUE..
  45. *          1 <= ILO <= max(1,IHI); IHI <= N.
  46. *
  47. *  H       (input/output) COMPLEX*16 array, dimension (LDH,N)
  48. *          On entry, the upper Hessenberg matrix H.
  49. *          On exit, if WANTT is .TRUE., H is upper triangular in rows
  50. *          and columns ILO:IHI, with any 2-by-2 diagonal blocks in
  51. *          standard form. If WANTT is .FALSE., the contents of H are
  52. *          unspecified on exit.
  53. *
  54. *  LDH     (input) INTEGER
  55. *          The leading dimension of the array H. LDH >= max(1,N).
  56. *
  57. *  W       (output) COMPLEX*16 array, dimension (N)
  58. *          The computed eigenvalues ILO to IHI are stored in the
  59. *          corresponding elements of W. If WANTT is .TRUE., the
  60. *          eigenvalues are stored in the same order as on the diagonal
  61. *          of the Schur form returned in H, with W(i) = H(i,i).
  62. *
  63. *  ILOZ    (input) INTEGER
  64. *  IHIZ    (input) INTEGER
  65. *          Specify the rows of Z to which transformations must be
  66. *          applied if WANTZ is .TRUE..
  67. *          1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
  68. *
  69. *  Z       (input/output) COMPLEX*16 array, dimension (LDZ,N)
  70. *          If WANTZ is .TRUE., on entry Z must contain the current
  71. *          matrix Z of transformations accumulated by ZHSEQR, and on
  72. *          exit Z has been updated; transformations are applied only to
  73. *          the submatrix Z(ILOZ:IHIZ,ILO:IHI).
  74. *          If WANTZ is .FALSE., Z is not referenced.
  75. *
  76. *  LDZ     (input) INTEGER
  77. *          The leading dimension of the array Z. LDZ >= max(1,N).
  78. *
  79. *  INFO    (output) INTEGER
  80. *          = 0: successful exit
  81. *          > 0: if INFO = i, ZLAHQR failed to compute all the
  82. *               eigenvalues ILO to IHI in a total of 30*(IHI-ILO+1)
  83. *               iterations; elements i+1:ihi of W contain those
  84. *               eigenvalues which have been successfully computed.
  85. *
  86. *  =====================================================================
  87. *
  88. *     .. Parameters ..
  89.       COMPLEX*16         ZERO, ONE
  90.       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ),
  91.      $                   ONE = ( 1.0D+0, 0.0D+0 ) )
  92.       DOUBLE PRECISION   RZERO, RONE, HALF
  93.       PARAMETER          ( RZERO = 0.0D+0, RONE = 1.0D+0,
  94.      $                   HALF = 0.5D+0 )
  95. *     ..
  96. *     .. Local Scalars ..
  97.       INTEGER            I, I1, I2, ITN, ITS, J, K, L, M, NH, NZ
  98.       DOUBLE PRECISION   H10, H21, OVFL, RTEMP, S, SMLNUM, T2, TST1,
  99.      $                   ULP, UNFL
  100.       COMPLEX*16         CDUM, H11, H11S, H22, SUM, T, T1, TEMP, U, V2,
  101.      $                   X, Y
  102. *     ..
  103. *     .. Local Arrays ..
  104.       DOUBLE PRECISION   RWORK( 1 )
  105.       COMPLEX*16         V( 2 )
  106. *     ..
  107. *     .. External Functions ..
  108.       DOUBLE PRECISION   DLAMCH, DLAPY2, ZLANHS
  109.       COMPLEX*16         ZLADIV
  110.       EXTERNAL           DLAMCH, DLAPY2, ZLANHS, ZLADIV
  111. *     ..
  112. *     .. External Subroutines ..
  113.       EXTERNAL           DLABAD, ZCOPY, ZLARFG, ZSCAL
  114. *     ..
  115. *     .. Intrinsic Functions ..
  116.       INTRINSIC          ABS, DBLE, DCONJG, DIMAG, MAX, MIN, SQRT
  117. *     ..
  118. *     .. Statement Functions ..
  119.       DOUBLE PRECISION   CABS1
  120. *     ..
  121. *     .. Statement Function definitions ..
  122.       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
  123. *     ..
  124. *     .. Executable Statements ..
  125. *
  126.       INFO = 0
  127. *
  128. *     Quick return if possible
  129. *
  130.       IF( N.EQ.0 )
  131.      $   RETURN
  132.       IF( ILO.EQ.IHI ) THEN
  133.          W( ILO ) = H( ILO, ILO )
  134.          RETURN
  135.       END IF
  136. *
  137.       NH = IHI - ILO + 1
  138.       NZ = IHIZ - ILOZ + 1
  139. *
  140. *     Set machine-dependent constants for the stopping criterion.
  141. *     If norm(H) <= sqrt(OVFL), overflow should not occur.
  142. *
  143.       UNFL = DLAMCH( 'Safe minimum' )
  144.       OVFL = RONE / UNFL
  145.       CALL DLABAD( UNFL, OVFL )
  146.       ULP = DLAMCH( 'Precision' )
  147.       SMLNUM = UNFL*( NH / ULP )
  148. *
  149. *     I1 and I2 are the indices of the first row and last column of H
  150. *     to which transformations must be applied. If eigenvalues only are
  151. *     being computed, I1 and I2 are set inside the main loop.
  152. *
  153.       IF( WANTT ) THEN
  154.          I1 = 1
  155.          I2 = N
  156.       END IF
  157. *
  158. *     ITN is the total number of QR iterations allowed.
  159. *
  160.       ITN = 30*NH
  161. *
  162. *     The main loop begins here. I is the loop index and decreases from
  163. *     IHI to ILO in steps of 1. Each iteration of the loop works
  164. *     with the active submatrix in rows and columns L to I.
  165. *     Eigenvalues I+1 to IHI have already converged. Either L = ILO, or
  166. *     H(L,L-1) is negligible so that the matrix splits.
  167. *
  168.       I = IHI
  169.    10 CONTINUE
  170.       IF( I.LT.ILO )
  171.      $   GO TO 130
  172. *
  173. *     Perform QR iterations on rows and columns ILO to I until a
  174. *     submatrix of order 1 splits off at the bottom because a
  175. *     subdiagonal element has become negligible.
  176. *
  177.       L = ILO
  178.       DO 110 ITS = 0, ITN
  179. *
  180. *        Look for a single small subdiagonal element.
  181. *
  182.          DO 20 K = I, L + 1, -1
  183.             TST1 = CABS1( H( K-1, K-1 ) ) + CABS1( H( K, K ) )
  184.             IF( TST1.EQ.RZERO )
  185.      $         TST1 = ZLANHS( '1', I-L+1, H( L, L ), LDH, RWORK )
  186.             IF( ABS( DBLE( H( K, K-1 ) ) ).LE.MAX( ULP*TST1, SMLNUM ) )
  187.      $         GO TO 30
  188.    20    CONTINUE
  189.    30    CONTINUE
  190.          L = K
  191.          IF( L.GT.ILO ) THEN
  192. *
  193. *           H(L,L-1) is negligible
  194. *
  195.             H( L, L-1 ) = ZERO
  196.          END IF
  197. *
  198. *        Exit from loop if a submatrix of order 1 has split off.
  199. *
  200.          IF( L.GE.I )
  201.      $      GO TO 120
  202. *
  203. *        Now the active submatrix is in rows and columns L to I. If
  204. *        eigenvalues only are being computed, only the active submatrix
  205. *        need be transformed.
  206. *
  207.          IF( .NOT.WANTT ) THEN
  208.             I1 = L
  209.             I2 = I
  210.          END IF
  211. *
  212.          IF( ITS.EQ.10 .OR. ITS.EQ.20 ) THEN
  213. *
  214. *           Exceptional shift.
  215. *
  216.             T = ABS( DBLE( H( I, I-1 ) ) ) +
  217.      $          ABS( DBLE( H( I-1, I-2 ) ) )
  218.          ELSE
  219. *
  220. *           Wilkinson's shift.
  221. *
  222.             T = H( I, I )
  223.             U = H( I-1, I )*DBLE( H( I, I-1 ) )
  224.             IF( U.NE.ZERO ) THEN
  225.                X = HALF*( H( I-1, I-1 )-T )
  226.                Y = SQRT( X*X+U )
  227.                IF( DBLE( X )*DBLE( Y )+DIMAG( X )*DIMAG( Y ).LT.RZERO )
  228.      $            Y = -Y
  229.                T = T - ZLADIV( U, ( X+Y ) )
  230.             END IF
  231.          END IF
  232. *
  233. *        Look for two consecutive small subdiagonal elements.
  234. *
  235.          DO 40 M = I - 1, L, -1
  236. *
  237. *           Determine the effect of starting the single-shift QR
  238. *           iteration at row M, and see if this would make H(M,M-1)
  239. *           negligible.
  240. *
  241.             H11 = H( M, M )
  242.             H22 = H( M+1, M+1 )
  243.             H11S = H11 - T
  244.             H21 = H( M+1, M )
  245.             S = CABS1( H11S ) + ABS( H21 )
  246.             H11S = H11S / S
  247.             H21 = H21 / S
  248.             V( 1 ) = H11S
  249.             V( 2 ) = H21
  250.             IF( M.EQ.L )
  251.      $         GO TO 50
  252.             H10 = H( M, M-1 )
  253.             TST1 = CABS1( H11S )*( CABS1( H11 )+CABS1( H22 ) )
  254.             IF( ABS( H10*H21 ).LE.ULP*TST1 )
  255.      $         GO TO 50
  256.    40    CONTINUE
  257.    50    CONTINUE
  258. *
  259. *        Single-shift QR step
  260. *
  261.          DO 100 K = M, I - 1
  262. *
  263. *           The first iteration of this loop determines a reflection G
  264. *           from the vector V and applies it from left and right to H,
  265. *           thus creating a nonzero bulge below the subdiagonal.
  266. *
  267. *           Each subsequent iteration determines a reflection G to
  268. *           restore the Hessenberg form in the (K-1)th column, and thus
  269. *           chases the bulge one step toward the bottom of the active
  270. *           submatrix.
  271. *
  272. *           V(2) is always real before the call to ZLARFG, and hence
  273. *           after the call T2 ( = T1*V(2) ) is also real.
  274. *
  275.             IF( K.GT.M )
  276.      $         CALL ZCOPY( 2, H( K, K-1 ), 1, V, 1 )
  277.             CALL ZLARFG( 2, V( 1 ), V( 2 ), 1, T1 )
  278.             IF( K.GT.M ) THEN
  279.                H( K, K-1 ) = V( 1 )
  280.                H( K+1, K-1 ) = ZERO
  281.             END IF
  282.             V2 = V( 2 )
  283.             T2 = DBLE( T1*V2 )
  284. *
  285. *           Apply G from the left to transform the rows of the matrix
  286. *           in columns K to I2.
  287. *
  288.             DO 60 J = K, I2
  289.                SUM = DCONJG( T1 )*H( K, J ) + T2*H( K+1, J )
  290.                H( K, J ) = H( K, J ) - SUM
  291.                H( K+1, J ) = H( K+1, J ) - SUM*V2
  292.    60       CONTINUE
  293. *
  294. *           Apply G from the right to transform the columns of the
  295. *           matrix in rows I1 to min(K+2,I).
  296. *
  297.             DO 70 J = I1, MIN( K+2, I )
  298.                SUM = T1*H( J, K ) + T2*H( J, K+1 )
  299.                H( J, K ) = H( J, K ) - SUM
  300.                H( J, K+1 ) = H( J, K+1 ) - SUM*DCONJG( V2 )
  301.    70       CONTINUE
  302. *
  303.             IF( WANTZ ) THEN
  304. *
  305. *              Accumulate transformations in the matrix Z
  306. *
  307.                DO 80 J = ILOZ, IHIZ
  308.                   SUM = T1*Z( J, K ) + T2*Z( J, K+1 )
  309.                   Z( J, K ) = Z( J, K ) - SUM
  310.                   Z( J, K+1 ) = Z( J, K+1 ) - SUM*DCONJG( V2 )
  311.    80          CONTINUE
  312.             END IF
  313. *
  314.             IF( K.EQ.M .AND. M.GT.L ) THEN
  315. *
  316. *              If the QR step was started at row M > L because two
  317. *              consecutive small subdiagonals were found, then extra
  318. *              scaling must be performed to ensure that H(M,M-1) remains
  319. *              real.
  320. *
  321.                TEMP = ONE - T1
  322.                TEMP = TEMP / DLAPY2( DBLE( TEMP ), DIMAG( TEMP ) )
  323.                H( M+1, M ) = H( M+1, M )*DCONJG( TEMP )
  324.                IF( M+2.LE.I )
  325.      $            H( M+2, M+1 ) = H( M+2, M+1 )*TEMP
  326.                DO 90 J = M, I
  327.                   IF( J.NE.M+1 ) THEN
  328.                      IF( I2.GT.J )
  329.      $                  CALL ZSCAL( I2-J, TEMP, H( J, J+1 ), LDH )
  330.                      CALL ZSCAL( J-I1, DCONJG( TEMP ), H( I1, J ), 1 )
  331.                      IF( WANTZ ) THEN
  332.                         CALL ZSCAL( NZ, DCONJG( TEMP ), Z( ILOZ, J ),
  333.      $                              1 )
  334.                      END IF
  335.                   END IF
  336.    90          CONTINUE
  337.             END IF
  338.   100    CONTINUE
  339. *
  340. *        Ensure that H(I,I-1) is real.
  341. *
  342.          TEMP = H( I, I-1 )
  343.          IF( DIMAG( TEMP ).NE.RZERO ) THEN
  344.             RTEMP = DLAPY2( DBLE( TEMP ), DIMAG( TEMP ) )
  345.             H( I, I-1 ) = RTEMP
  346.             TEMP = TEMP / RTEMP
  347.             IF( I2.GT.I )
  348.      $         CALL ZSCAL( I2-I, DCONJG( TEMP ), H( I, I+1 ), LDH )
  349.             CALL ZSCAL( I-I1, TEMP, H( I1, I ), 1 )
  350.             IF( WANTZ ) THEN
  351.                CALL ZSCAL( NZ, TEMP, Z( ILOZ, I ), 1 )
  352.             END IF
  353.          END IF
  354. *
  355.   110 CONTINUE
  356. *
  357. *     Failure to converge in remaining number of iterations
  358. *
  359.       INFO = I
  360.       RETURN
  361. *
  362.   120 CONTINUE
  363. *
  364. *     H(I,I-1) is negligible: one eigenvalue has converged.
  365. *
  366.       W( I ) = H( I, I )
  367. *
  368. *     Decrement number of remaining iterations, and return to start of
  369. *     the main loop with new value of I.
  370. *
  371.       ITN = ITN - ITS
  372.       I = L - 1
  373.       GO TO 10
  374. *
  375.   130 CONTINUE
  376.       RETURN
  377. *
  378. *     End of ZLAHQR
  379. *
  380.       END
  381.